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ABSTRACT 

The  analysis  of  stresses  induced  by  contact  between  two  bodies  is  inherently  difficult 
because  the  size  of  the  contact  zone  is  unknown  and  constantly  changing  throughout 
loading.  To  overcome  these  difficulties,  two  approximation  methods  have  been  devel- 
oped to  determine  the  magnitude  of  contact  stresses  using  the  Rayleigh-Ritz  method  and 
the  finite  element  method.  Numerical  optimization  methods  are  employed  to  solve  the 
contact  problem.  The  solution  techniques  are  compared  to  known  analytical  solutions 
and  shown  to  yield  accurate  results.  An  application  of  this  approach  to  solving  the 
contact  problem  is  illustrated  by  examining  the  response  of  a  clamped  sandwich  com- 
posite beam  to  low  velocity  impact.  It  was  found  that  the  maximum  shear  stress  is  in- 
sensitive to  lamina  thickness,  however  an  increase  in  the  contact  layer  thickness  resulted 
in  a  reduction  in  interfacial  shear  stress.  In  addition,  it  was  noted  that  a  nonlinear 
bending  stress  distribution  in  the  contact  layer  intensified  as  the  thickness  of  this  layer 
increased.  This  phenomenon  was  found  to  be  localized  to  the  region  of  contact.  Finally, 
it  was  found  that  the  compressive  transverse  normal  stresses  increased  as  the  thickness 
of  the  contact  lamina  increased. 
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I.     INTRODUCTION 

A.     MOTIVATION 

Contact  stresses  occur  when  two  bodies  exert  forces  over  limited  contact  regions. 
Examples  of  contact  stress  include  meshing  gear  teeth,  cam  shaft  and  pushrod  contact, 
rollers  in  plate  forming  operations,  roller  and  ball  bearings  in  contact  with  races,  shaft 
and  journal  bearing  contact,  and  plate-pin  connections. 

Contact  zones  can  be  point,  line,  or  surface  in  geometry.  Because  of  the  limited 
contact  zone,  the  local  stresses  can  be  sufficiently  high  to  be  of  major  concern  to  the 
designer.  Consequently,  a  thorough  understanding  of  this  phenomenon  is  essential.  The 
first  successful  analytical  solution  to  the  contact  problem  between  two  spheres  was 
solved  in  the  late  19th  century  by  Hertz.  His  solution  can  only  be  applied  to  simple 
geometries  such  as  spheres,  cylinders,  and  flat  plates.  Because  of  these  limitations,  al- 
ternate solution  techniques  were  needed  to  accommodate  more  complex  geometries  and 
boundary  conditions. 

Unfortunately,  the  contact  problem  is  difficult  to  study.  The  most  significant  diffi- 
culty is  that  the  size  of  the  contact  zone  is  unknown  and  constantly  changing  throughout 
loading.  Solutions  are  obtained  by  an  iterative  process.  Consequently,  the  problem  is 
highly  nonlinear  in  its  behavior.  It  is  because  of  these  difficulties  that  approximation 
and  numerical  methods  are  preferred  in  the  solution  of  the  contact  problem.  There  are 
two  general  approaches  used  in  the  approximate  techniques  to  solve  the  contact  prob- 
lem. One  class  applies  specific  iterative  procedures  to  solve  a  nonlinear  system  of 
equations  that  represent  the  contact  state.  For  example,  the  contact  condition  can  be 
simulated  by  the  introduction  of  additional  coupling  terms  into  the  system  of  equations. 
The  other  class  constructs  a  functional  that  includes  the  contact  body  constraint.  The 
functional  is  then  minimized  using  specific  numerical  programming  techniques.  Al- 
though these  two  classes  use  different  procedures,  there  are  several  methods  to  incor- 
porate the  contact  condition  into  the  problem  formulation  that  are  common  to  both 
classes.  Examples  of  these  include  the  penalty  method  and  the  augmented  Lagrange 
multiplier  method.  These  methods  specify  the  manner  in  which  the  contact  boundary 
conditions  are  treated. 

The  objective  of  this  study  has  two  parts.  First,  two  approximate  solution  tech- 
niques will  be  developed  to  obtain  the  solution  of  the  general  contact  problem.    These 


techniques  belong  to  the  general  class  of  methods  that  calculate  a  specific  functional  and 
then  applies  numerical  programming  methods  to  solve  the  contact  problem.  Second, 
these  solution  techniques  will  be  verified  by  comparison  with  known  analytical  solutions. 
With  the  solution  techniques  verified,  the  methods  will  be  applied  to  study  an  actual 
contact  problem. 

In  order  to  accomplish  the  first  objective,  two  models  will  be  developed.  The  first 
method  utilizes  the  Rayleigh-Ritz  method  to  solve  the  contact  problem.  An  assumed 
deformation  field  that  satisfies  the  boundary  conditions  is  carefully  selected.  Theory  of 
elasticity  relationships  are  applied  to  obtain  an  expression  for  the  system's  strain  energy 
enabling  calculation  of  the  total  potential.  The  minimization  of  the  total  potential  en- 
ergy enables  the  calculation  of  the  contact  stresses  at  any  point  in  the  body.  The  second 
model  developed  uses  finite  element  analysis  to  accomplish  the  same  objectives.  The 
numerical  minimization  technique  used  in  both  cases  is  the  augmented  Lagrange  multi- 
plier method. 

To  accomplish  the  second  objective,  both  methods  will  be  verified  by  comparison 
with  the  Hertz  solution  of  an  infinitely  long  cylinder  in  contact  with  a  flat  plane.  With 
verification  completed,  the  contact  stresses  resulting  from  low  velocity  impact  between 
objects  and  composite  sandwich  materials  will  be  studied. 

B.     LITERATURE  SURVEY 

As  discussed  previously,  there  are  two  general  approaches  to  solving  the  contact 
problem.  One  approach  uses  special  procedures  to  solve  a  nonlinear  system  of 
equations.  The  other  creates  a  functional  that  includes  the  contact  boundary  con- 
straints. Although  different,  both  approaches  use  similar  mathematical  methods  to  in- 
corporate the  contact  condition  into  the  problem  formulation.  Two  of  these  methods 
are  the  Lagrange  method  and  the  penalty  method.  A  detailed  explanation  of  the  method 
of  employing  the  Lagrange  method  in  nonlinear  finite  element  analysis  was  discussed  by 
Pian  [Ref.  1].  Alternately,  the  penalty  finite  element  method  was  used  by  Cheng  [Ref. 
2]  to  solve  the  multibody  contact  problem.  The  fundamental  concept  of  the  latter  ap- 
proach is  the  transformation  of  a  constrained  problem  into  an  unconstrained  one.  This 
is  done  by  penalizing  (i.e.,  increasing)  the  objective  function  for  constraint  violations. 
Although  both  methods  are  effective,  there  is  substantial  discussion  on  the  limitations 
of  both  methods.  Nour-Omid  [Ref.  3]  described  the  positive  and  negative  aspects  of 
both  methods.  The  Lagrange  method  has  been  shown  to  be  the  more  accurate  method. 
However,  its  usage  requires  the  introduction  of  additional  unknowns  thus  increasing  the 


II.     FORMULATION  OF  THE  CONTACT  PROBLEM 

A.     PRINCIPLE  OF  MINIMUM  TOTAL  POTENTIAL  ENERGY 

As  discussed  in  the  introduction,  the  limited  utility  of  the  analytical  solutions  ne- 
cessitated the  development  of  solution  techniques  capable  of  handling  the  nonlinear  be- 
havior of  the  contact  problem  with  complicated  geometry  and  complex  boundary 
conditions.  This  study  intends  to  develop  two  numerical  procedures  to  solve  the  contact 
problem.  In  short,  the  procedures  will  use  different  methods  to  obtain  a  functional,  the 
system's  total  potential  energy,  and  then  use  similar  methods  to  obtain  the  equilibrium 
condition.  Determination  of  the  equilibrium  position  is  made  by  application  of  the 
principle  of  minimum  potential  energy.  With  equilibrium  established,  contact  stresses 
can  be  quantified.  In  order  to  understand  the  details  of  this  approach,  the  principle  of 
virtual  work  and  the  principle  of  minimum  potential  energy  must  be  discussed. 

Given  a  body  in  equilibrium,  it  is  desired  to  describe  the  response  of  that  body  to 
infinitesimal  displacements  resulting  from  a  system  of  forces.  If  each  particle  in  the  body 
is  described  by  some  generalized  coordinates,  then  the  work  resulting  from  these 
infinitesimal  displacements  is  simply  the  product  of  the  generalized  forces  acting  on  each 
particle  and  the  particle's  displacement.  However,  if  the  particle  is  in  equilibrium,  the 
work  must  be  zero  since  the  summation  of  forces  in  the  x,  y,  and  z  directions  is  zero. 
The  infinitesimal  displacements  and  work  in  this  example  are  referred  to  as  virtual  in 
nature.    The  fact  that  this  work  vanishes  is  referred  to  as  the  principle  of  virtual  work. 

The  virtual  work  discussed  thus  far  can  be  subcategorized  as  virtual  strain  energy 
and  virtual  work  done  by  external  forces.  From  the  definition  of  strain  energy,  virtual 
strain  energy,  SU,  that  results  from  virtual  displacements  can  be  calculated.  Since  this 
energy  is  viewed  as  energy  against  the  bonds  between  elements,  SU  is  a  negative  quan- 
tity. The  work  done  by  external  forces  is  designated  S  XV  and  is  simply  the  summation 
of  the  product  of  the  external  forces  and  the  displacements  of  the  generalized  coordi- 
nates. 

Since  the  principle  of  virtual  work  states  that  the  work  done  as  a  result  of  virtual 
displacements  is  zero, 

SlV-SU  =  0 


Alternately,  this  can  be  expressed  as, 

n  =d(U-  ll')  =  0 
where  n  represents  the  system's  total  potential. 

The  above  equation  illustrates  the  condition  of  minimum  total  potential  of  a  system. 
This  is  the  foundation  of  the  principle  of  minimum  potential  energy.  [Ref.  12:  pp. 
330-331]  This  principle  states  that  in  a  condition  of  stable  equilibrium,  the  system's 
total  potential  is  stationary.  Hence,  determination  of  a  system's  total  potential  energy 
and  the  minimization  of  that  quantity  will  enable  the  calculation  of  the  equilibrium  po- 
sition.   This  is  the  basis  for  the  numerical  techniques  developed  in  this  study. 

B.     CONTACT  PROBLEM  DESCRIPTION 

The  objective  of  the  solution  techniques  to  be  developed  is  to  obtain  a  means  of 
calculating  the  total  potential  energy  and  minimize  it  to  determine  the  equilibrium  con- 
dition. To  accomplish  this,  two  different  models  will  be  created  to  first  solve  a  simple 
isotropic  case,  the  solution  of  which  is  known.  In  this  manner,  our  models  can  be  vali- 
dated for  use  in  more  complicated  arrangements. 

Consider  contact  between  a  cylinder  and  an  infinite  plane  as  shown  in  Figure  1.  A 
cross  section  of  the  contact  zone  is  shown  in  Figure  2.  Let  the  width  of  the  contact  zone 
equal  a  distance  of  2a.  An  analytical  solution  of  this  problem  is  available  as  a  result  of 
the  work  done  by  Hertz.  In  developing  the  solution  methods,  only  a  very  limited  region 
adjacent  to  the  contact  zone  will  be  examined.  The  reason  for  this  is  that  the  contact 
phenomenon  is  a  very  local  one.  Figure  3  represents  the  analytical  solution  of  the 
normal  stresses  resulting  from  this  contact  problem.  [Ref.  13]  As  shown,  the  stresses 
in  the  foundation  diminish  very  rapidly.  The  diagram  shows  that  az  becomes  negligible 
at  depths  less  than  one  half-contact  zone  (a)  and  ay  diminishes  significantly  in  less  than 
3  half  zones.  Because  of  this,  it  is  reasonable  to  assume  that  displacements  beyond  a 
very  limited  region  are  negligible  in  the  strain  energy  calculation  of  the  contact  problem. 

A  number  of  simplifying  assumptions  will  be  made  for  this  study: 

1.  As  discussed  above,  displacements  beyond  a  limited  region  are  negligible  in  strain 
energy  calculations. 

2.  The  foundation  is  an  elastic  isotropic  material.   The  cylinder  (roller)  is  rigid. 


system's  degrees-of-freedom  and  computational  time.  The  penalty  method  does  not  re- 
quire additional  computational  time,  but  it  has  been  shown  to  result  in  less  accurate 
solutions  in  satisfying  the  contact  boundary  conditions.  Since  the  contact  boundary 
conditions  are  exactly  satisfied  only  when  the  penalty  term  goes  to  infinity,  the  correct 
choice  of  the  penalty  parameter  is  the  key  to  an  acceptable  solution.  Guerra  [Ref.  4] 
supported  these  claims. 

Because  of  the  limitations  o[  the  Lagrange  and  penalty  methods,  Bishoff  [Ref.  5] 
advocated  the  use  of  the  augmented  Lagrange  multiplier  method  to  solve  finite  element 
contact  problems.  This  method  is  favorable  since  it  avoids  the  limitations  of  both 
methods  discussed  above.  The  augmented  Lagrange  multiplier  method  is  similar  to  the 
penalty  method  in  that  the  objective  function  is  penalized  for  constraint  violations. 
However,  an  additional  multiplier  term  is  added  so  the  optimum  can  be  achieved  by  a 
combination  of  both  terms.  As  stated  by  Vanderplaats  [Ref.  6:  p.  141],  the  advantage 
of  this  method  is  that  the  penalty  term  is  not  required  to  grow  to  infinity  to  achieve  exact 
constraint  satisfaction. 

The  augmented  Lagrange  multiplier  method  can  be  used  in  numerical  programming 
techniques  to  find  the  optimum  of  any  functional.  Pierre  and  Lowe  [Ref.  7]  provide  a 
detailed  analysis  of  the  programming  techniques  necessary  in  applying  this  method. 
Additionally,  Vanderplaats  [Ref.  6:  pp.  140-147]  provided  an  excellent  discussion  on  the 
practical  usage  of  this  method  in  numerical  techniques.  Rothert  et  al.  [Ref.  8]  used  a 
numerical  programming  code  based  on  this  method  to  solve  a  nonlinear  contact  prob- 
lem. In  this  study,  an  existing  numerical  optimization  routine  utilizing  the  augmented 
Lagrange  multiplier  method  will  be  used  as  an  integral  part  of  two  solution  methods 
developed  to  solve  the  contact  problem.  These  methods  will  use  different  techniques  to 
obtain  the  same  functional.  The  augmented  Lagrange  multiplier  method  will  then  be 
used  in  a  similar  fashion  to  solve  each  problem.  In  each  case,  the  optimization  routine 
will  be  used  to  determine  a  set  of  design  variables  that  describes  the  contact  state. 

Following  the  development  and  verification  of  the  numerical  procedures,  this  study 
will  investigate  the  response  of  composite  sandwich  materials  to  low  velocity  impact. 
One  of  the  common  failure  mode  of  low  velocity  impact  is  delamination.  Joshi  and  Sun 
[Ref.  9]  studied  the  impact  response  of  a  three  layer  cross-ply  graphite  epoxy  laminate. 
A  correlation  was  obtained  between  delamination  cracks  initiated  experimentally  and 
maximum  shear  stress  points  determined  numerically.  Sun  and  Rechak  [Ref.  10]  fol- 
lowed up  these  findings  and  found  that  the  introduction  of  adhesive  layers  between 
laminae  reduced  the  shear  stress  distribution  thus  reducing  delamination.    Choi,  Wang, 


and  Chang  [Ref.  11]  studied  the  effects  of  laminae  orientation,  ply  thickness,  and 
stacking  sequence  on  impact  damage  of  graphite  epoxy  composites.  It  was  determined 
that  stacking  sequence  affects  impact  damage  more  than  laminae  thickness  variations. 
Much  of  the  previous  work  has  focused  on  the  behavior  of  the  graphite  epoxy  laminate. 
However,  there  is  currently  interest  in  the  development  of  turbine  blades  constructed  of 
sandwich  composites.  It  is  therefore  beneficial  to  investigate  the  response  of  composite 
sandwich  materials  to  low  velocity  impact. 


Figure   1.      Roller- foundation  assembly 


Figure  2.      Contact  zone  cross  section 


3.  Deformations  normal  to  the  cross-section  are  negligible,  hence  a  condition  of  plane 
strain  exists. 

4.  The  roller  is  subjected  to  a  vertical  distributed  load. 

5.  The  roller-foundation  contact  is  frictionless. 

An  important  restriction  on  the  minimization  problem  will  be  that  one  body  will  be 
prohibited  from  penetrating  into  the  other  body.  This  may  seem  like  an  obvious  re- 
striction, however  a  method  of  mathematically  stating  this  restriction  must  be  discussed. 
Figure  4  shows  the  deformed  and  undeformed  contact  zone.  Let  <5  represent  the  de- 
flection of  the  roller  due  to  an  external  force  and  v(jCjy)  represent  the  vertical  deflection 
of  the  foundation  at  any  point  (xj,-).  At  the  point  of  Contact  A,  the  condition  of  no 
interference  can  be  expressed  as. 

v(0.0)  >  c) 

At  point  B,  this  condition  can  be  stated  as. 

\{r  sin  0,0)  >d  -  r{\  -  cos  6) 

The  latter  condition  can  be  specified  at  as  many  points  as  necessary  to  define  this  re- 
striction. 

C.     NUMERICAL  OPTIMIZATION 

1.     Optimization  Fundamentals 

Before  developing  the  models  to  be  used  in  this  study,  one  final  area  must  be 
discussed.  Once  the  total  potential  energy  has  been  calculated,  a  method  of  minimizing 
it  to  find  the  equilibrium  position  must  be  employed.  A  study  of  the  numerical  opti- 
mization technique  to  be  used  is  required. 

The  technique  being  used  in  this  study  relies  heavily  on  the  methods  of  design 
optimization.  Design  optimization  is  the  utilization  of  mathematical  techniques  to 
minimize  or  maximize  a  particular  value  to  obtain  an  optimum  solution.  The  method 
is  ideally  suited  for  design.  A  given  design  task  may  have  an  infinite  number  of  sol- 
utions. However,  finding  the  best  solution  is  a  matter  of  the  designer's  experience  and 
intuition.  In  the  absence  of  significant  experience  in  a  particular  field,  finding  this  sol- 
ution may  reduce  to  examining  a  range  of  possible  solutions  by  trial  and  error.  Opti- 
mization routines  can  be  utilized  to  find  this  solution  mathematically. 

Optimization  problems  can  be  constrained  or  unconstrained.  For  example,  it 
may  be  desired  to  determine  the  minimum  of  the  parabola,  F{x)  =  (x  -  5)2  +  2.   As  seen 
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Figure  3.      Analytical  solutions  for  ay  and  <x, 


in  Figure  5,  the  minimum  is  clearly  identified  at  point  A.  This  is  an  example  of  an  un- 
constrained problem.   A  constrained  counterpart  of  this  problem  is: 

Minimize:     F{x) 
Subject  to:  ^(.v)  >  .5jc  +  4 

As  seen  in  Figure  6,  the  minimum  of  the  constrained  problem  is  at  point  B.  This  is  a 
simple  illustration  of  constrained  minimization.  The  contact  problem  is  a  constrained 
minimization  problem. 


Figure  4.      Deformed  contact  zone 

The  value  to  be  minimized  or  maximized  is  referred  to  as  the  objective  function. 
The  parameters  to  be  determined  are  referred  to  as  design  variables.  The  optimization 
process  is  an  iterative  one.  The  objective  function  is  evaluated.  The  design  variables 
are  varied  thus  obtaining  a  new  objective  function.  If  the  difference  is  within  a  certain 
tolerance  or  meets  a  certain  convergence  criteria,  the  optimum  has  been  obtained.  If  it 
is  out  of  tolerance,  the  cycle  repeats. 
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Figure  5.      Unconstrained  minimization 
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2.     Augmented  Lagrange  Multiplier  Method 

The  contact  problem  belongs  in  the  class  of  constrained  problems.  There  are 
several  techniques  of  solving  constrained  problems.  The  technique  used  in  this  study 
belongs  to  a  class  of  solution  techniques  known  as  sequential  unconstrained  minimiza- 
tion techniques  (SUMT).  This  class  of  techniques  is  designed  for  the  general  nonlinear 
problem.  The  fundamental  concept  behind  this  approach  is  that  a  constrained  problem 
is  transformed  into  an  unconstrained  problem  and  the  objective  function  is  minimized 
using  an  unconstrained  minimization  technique.  This  transformation  is  accomplished 
by  assessing  a  penalty  to  the  objective  function  for  constraint  violations.  For  example, 
if  the  design  variables  are  varied  in  such  a  way  as  to  enter  the  region  where  a  constraint 
is  violated  (i.e.,  the  infeasible  region),  the  objective  function  would  be  assessed  a  penalty 
(i.e.,  increased).  Thus,  in  order  to  minimize  the  objective  function,  the  design  variables 
remain  within  the  region  where  no  constraints  are  violated  (i.e..  the  feasible  region). 
[Ref.    6:  pp.  121-123] 

There  are  a  number  of  methods  within  this  class  of  techniques.  They  essentially 
differ  in  the  way  in  which  penalties  are  assessed.  The  technique  used  in  this  study  is  the 
augmented  Lagrange  multiplier  (ALM)  method. 

Given  the  constrained  inequality  optimization  problem; 

Minimize:     F{X) 

Subject  to:  gt(X)  <  0,  /  =  1 .2 n 

The  Augmented  Lagrangian  is  defined  as, 

m 

A(X,  k,p)  =  F[X)  +  £{/M  +  s,2]  +  p[g/LX)  +  sff} 

where, 

X=  vector  containing  the  design  variables 

).t  =  Lagrange  multipliers 

p  =penalty  parameter 

5,  =slack  variables  which  convert  inequality  constraints  to  equality  con- 
straints 
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The  first  two  terms  of  A(X,  X,p)  represent  the  Lagrangian.  From  the  method  of 
Lagrange  multipliers,  it  is  known  that  minimization  of  the  Lagrangian  represents  opti- 
mality.  Like  simple  problems  where  X  is  simply  an  additional  unknown  to  obtain,  in  the 
AL.M  method  X  is  unknown.  Hence,  a  mathematical  routine  based  on  constraint  values 
is  used  to  select  and  modify  /  for  successive  iterations.  [Ref.  6:  pp.  140-147]  Initial 
selection  of  this  term  can  have  a  significant  impact  of  the  problem's  convergence. 

As  discussed  above,  a  penalty  is  assessed  to  the  objective  function  for  constraint 
violations.  This  feature  is  apparent  by  examining  the  last  term.  A  constraint  violation 
results  in  a  positive  value  for  g(X)  thus  resulting  in  an  increase  in  A.  The  value  of/?  is 
a  scaling  term  which  is  sequentially  increased  throughout  optimization.  This  ensures 
that  there  is  a  balance  between  convergence  and  numerical  conditioning.  I  f  /?  remains 
small,  convergence  may  occur  with  major  constraint  violations.  If  p  remains  large, 
constraints  will  be  satisfied  at  the  expense  of  an  ill-conditioned  problem  [Ref.  7:  p. 
169].  As  with  the  Lagrange  term,  selection  of  this  term  has  a  significant  impact  on  the 
outcome  of  the  problem. 

3.     Optimizer  and  One-Dimensional  Search 

Thus  far.  a  procedure  has  been  defined  which  has  transformed  a  constrained 
minimization  problem  into  an  unconstrained  one.  This  level  of  the  optimization  process 
is  referred  to  as  the  optimization  strategy.  The  formulation  of  the  modified  objective 
function  via  the  augmented  Lagrange  multiplier  method  represents  a  key  portion  of  the 
optimization  process.  However,  there  are  additional  parts  of  this  process  that  require 
comment. 

With  the  modified  objective  function  and  a  procedure  for  assessing  penalties  in 
place,  a  procedure  for  minimizing  the  objective  function  must  be  defined.  This  portion 
of  the  process  is  carried  out  by  the  'optimizer.'  The  optimizer's  function  is  to  system- 
atically alter  the  design  variables  in  a  manner  that  reduces  the  objective  function  rapidly. 
If  X  represents  a  vector  containing  the  design  variables,  the  following  process  is  used  by 
the  optimizer  to  alter  X 


where, 


Xl  —  <*(/_!)  +  *<->< 


/  =  current  iteration  number 

5  =  search  vector 

k  =  scalar  representing  distance  traveled  in  direction  S 
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In  general,  two  processes  must  be  accomplished  to  find  the  optimum.  First,  the  search 
direction  S  must  be  determined  by  a  systematic  process.  An  example  of  this  phase  is  the 
steepest  descent  method  where  the  direction  of  steepest  gradient  is  chosen.  Second,  the 
scalar  k  must  be  determined  such  that  the  objective  function  is  minimized  as  much  as 
possible  in  the  search  direction  of  the  current  iteration.  The  latter  phase  is  referred  to 
as  one-dimensional  search.  [Ref.  6:  pp.  10-12] 

The  optimizer  used  in  this  study  is  the  variable  metric  method.  Due  to  the 
complexity  of  this  method,  a  discussion  of  their  formulation  is  omitted.  A  detailed  ac- 
count is  available  in  Vanderplaats  [Ref.  6:  pp.  92-93].  The  one-dimensional  search 
routine  used  in  this  study  is  the  golden  section  method  with  polynomial  interpolation. 
The  one-dimensional  search  portion  of  the  process  merely  represents  a  systematic  and 
efficient  method  of  finding  the  minimum  in  the  chosen  search  direction.  A  detailed  ac- 
count of  this  approach  is  again  available  in  Vanderplaats  [Ref.  6:  pp.  26-49]. 
4.     Convergence 

The  final  point  to  be  discussed  relevant  to  optimization  fundamentals  is  that  of 
convergence.  Convergence  criteria  are  utilized  to  identify  the  optimum  solution  and 
terminate  calculations.  There  are  a  number  of  convergence  criteria  that  can  be  utilized. 
The  most  obvious  is  absolute  convergence  where  the  objective  functions  from  two  suc- 
cessive iterations  are  compared.  If  the  difference  between  the  values  is  within  some 
prescribed  limit,  optimization  is  terminated.  A  second  method  signifying  optimality  for 
unconstrained  problems  is  calculation  of  the  gradient  with  respect  to  the  design  variable 
vector  A'.  If  this  value  is  approximately  zero,  optimality  has  been  achieved.  This  method 
is  called  the  Kuhn-Tucker  conditions  for  unconstrained  minimization.  Kuhn-Tucker 
conditions  are  more  involved  for  constrained  problems.  [Ref.  6:  pp.    100-101] 

The  Automated  Design  Synthesis  (ADS)  System  used  in  this  study  uses  both 
these  termination  criteria  as  well  as  relative  convergence.  Relative  convergence  is  similar 
to  absolute  convergence  except  normalized  versions  of  the  difference  between  successive 
iterations  is  calculated.  Again,  if  a  specified  tolerance  is  achieved,  optimization  is  ter- 
minated. 
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III.     APPROXIMATE  SOLUTION  TECHNIQUES 

A.     RAYLEIGH-RITZ  APPROACH 

1.     Background 

The  Rayleigh-Ritz  method  is  a  method  of  utilizing  the  theory  of  minimum  po 
tential  energy  to  solve  a  given  problem.  The  fundamental  concept  behind  this  method 
is  that  a  trial  function  that  represents  the  deformation  field  is  chosen  in  terms  oC  un- 
known constants.  Second,  the  system's  total  potential  energy  is  calculated  in  terms  ofl 
the  trial  function.  Since  the  total  potential  is  a  minimum  at  equilibrium,  minimization 
enables  determination  of  the  unknown  constants.  The  total  potential  is  minimized  by 
differentiating  with  respect  to  the  unknown  constants  and  equating  to  zero.  The  result 
is   n'  equations  and  'n'  unknowns,  the  trial  function  constants.  [Ref.  12:  pp.  335-336  ] 

The  only  requirement  of  the  Rayleigh-Ritz  method  is  that  the  trial  function  is 
kinematically  admissible.  A  kinematically  admissible  solution  is  one  that  satisfies  the 
geometric  boundary  conditions  of  the  system  (i.e.,  deflection  and  slope).  Other  require- 
ments need  not  be  satisfied.  As  an  example,  consider  a  simply  supported  beam  of  length 
L  with  the  origin  at  the  left  end  of  the  beam.  A  kinematically  admissible  solution  to 
describe  the  beam's  one  dimensional  deformation  from  vertical  loading  in  the  y  direction 
is, 

y(x)  =  an  sin(  ~-  ) 
where  an  represents  the  coefficient  to  be  determined. 

Deflection  boundary  conditions  at  x  =  0  and  jc  =  L  have  been  satisfied.  Obviously,  an 
increased  number  of  terms  in  the  trial  function  will  yield  a  far  more  accurate  solution. 
A  Fourier  sine  series  would  be  a  reasonable  selection  in  this  case. 

Although  the  Rayleigh-Ritz  method  does  not  stipulate  numerous  requirements 
on  the  trial  function,  sensible  choices  of  trial  functions  will  increase  solution  accuracy 
significantly.  For  example,  consider  the  same  simply  supported  beam  with  the  origin 
now  at  the  center.   A  kinematically  admissible  function  is, 

/     \  •     /    27T.V    , 

y[x)  =  an  sin(  — —  ). 
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Naturally  however,  due  to  the  placement  of  the  origin  in  this  problem,  an  even  function 
is  a  much  for  sensible  selection  for  a  trial  function.  A  more  appropriate  selection  would 
be, 

y(x)  =  an  cos(  -^-  ). 

The  latter  point  concerning  sensible  selection  of  the  trial  function  will  be  discussed  in 
detail  throughout  this  study. 

2.  Application  of  the  Rayleigh-Ritz  Method  to  the  Contact  Problem. 

A  brief  overview  of  the  method  to  be  developed  is  in  order.  Application  of  the 
Rayleigh-Ritz  method  necessitates  the  selection  of  the  appropriate  trial  function  in  terms 
of  unknown  coefficients.  A  discussion  of  the  physical  nature  of  this  problem  as  well  as 
the  desired  features  of  the  trial  solution  is  required. 

The  theory  of  elasticity  relationships  will  be  applied  using  the  trial  function  to 
obtain  the  system's  strain  energy  in  terms  of  unknown  coefficients.  The  system's  total 
potential  energy  will  then  be  minimized  utilizing  the  optimization  techniques  discussed 
in  Chapter  II.  The  design  variables  are  the  unknown  trial  function  coefficients.  With 
the  coefficients  determined,  the  displacement  is  known  for  all  points  enabling  the  stress 
to  be  determined  throughout  the  body.  Figure  7  is  a  flow  chart  of  the  procedure  to  be 
utilized.  The  post-processing  procedure  shown  is  simply  the  calculation  of  the  stresses 
using  the  now  determined  coefficients. 

3.  Trial  Function  Selection 

In  order  to  choose  an  appropriate  trial  function,  it  is  necessary  to  have  an 
understanding  of  the  physical  phenomena  to  be  modeled.  Consider  the  roller-foundation 
system  shown  in  Figure  1.  Two  trial  functions  are  needed  to  model  horizontal  and  ver- 
tical deformation.  There  are  a  number  of  important  characteristics  that  should  be  in- 
herent within  the  trial  functions.   These  are  outlined  below: 

1.  As  seen  in  Figure  3,  at  and  ay  are  equal  and  compressive  at  the  point  of  contact. 
Additionally,  at  decays  more  rapidly  than  ay  as  distance  from  the  contact  point  in- 
creases. 

2.  Taking  the  origin  at  the  point  of  contact  as  shown  in  Figure  2,  the  u-deformation 
takes  the  form  of  an  odd  function. 

3.  The  v-deformation  takes  the  form  of  an  even  function  (i.e.  symmetric  about  the 
origin  and  greatest  at  the  origin). 

4.  To  satisfy  the  Rayleigh-Ritz  method  requirements,  the  boundary  conditions  must 
be  satisfied.  In  this  case,  this  requires  u  and  v  deformation  to  be  zero  at  the  far 
boundaries. 
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The  characteristics  outlined  in  item  4  satisfy  the  requirements  of  this  method. 
However,  since  this  study  will  be  a  stress  analysis,  it  is  also  desired  that  the  stresses 
also  reflect  the  physical  phenomena.  Since  as  and  ay  are  functions  of  ct  and  ct.  c, 
and  cy  should  also  exhibit  certain  characteristics.  Referring  to  Figure  8  for  dimen- 
sions and  the  coordinate  system,  it  is  desired  that  cx  and  cy  equal  zero  at  x=  L  2  and 
y=  H.  This  will  ensure  that  stresses  are  zero  at  the  boundaries.  To  satisfy  the  re- 
quirements of  item  1  above,  c,  and  cy  should  be  equal  and  negative  at  the  point  of 
contact  and  decrease  in  magnitude  as  the  distance  from  the  point  of  contact  in- 
creases. 


Figure  8.      Contact  zone 

With  the  above  guidelines  in  mind,  the  trial  function  can  be  selected.    The  trial 
functions  chosen  for  this  study  are  composed  of  a  series  of  terms  of  the  general  form: 


u(x,y)  =  an{H-y)a(  —  -x) 


(3- 


v(xiy)  =  bn(H-y)c(j--x)c 


(3.2) 


These  expressions  were  carefully  chosen  and  represent  a  compromise  due  to  the  diffi- 
culties of  satisfying  all  boundary  conditions  with  the  physical  phenomena  of  this  prob- 
lem. 
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From  the  above  expressions,  it  is  immediately  obvious  that  the  geometric 
boundary  conditions  at  x=L2  and  y=H  have  been  satisfied.  This  satisfies  the  re- 
quirements of  the  Rayleigh-Ritz  method.  In  addition,  there  are  a  number  of  important 
characteristics  that  illustrate  the  advantage  of  this  selection: 

1.  ix  and  £y  are  negative  thus  simulating  a  compressive  environment  in  the  vicinity  of 
the  point  of  contact.  The  importance  of  this  is  obvious.  If  normal  strains  were  not 
negative,  the  resulting  requirement  would  be  for  the  coefficient  to  be  less  than  zero 
to  simulate  compression.  It  is  obvious  that  this  would  result  in  deformations  op- 
posite to  that  which  was  desired  by  the  choice  of  the  trial  function. 

2.  This  selection  for  deformation  fields  has  the  important  characteristic  of  decreasing, 
deformation  as  we  move  away  from  the  point  of  contact.  Also  note  that  defor- 
mation is  maximum  at  the  point  of  contact. 

3.  The  exponents  a,b.c.  and  d  can  be  varied  to  simulate  subsurface  stress  fields.  For 
example,  if  the  analytical  solution  indicates  a  large  y  gradient  for  a,  at  x  =  0,  the: 
objective  would  be  to  increase  the  rate  at  which  c,  decreases  as  the  distance  from 
the  point  of  contact  increases.  This  would  be  easily  simulated  by  raising  the  value: 
of  a.  If  this  change  had  a  detrimental  effect  on  the  behavior  of  av,  the  exponents 
of  the  vertical  deformation  could  be  varied  to  restore  the  solution. 

It  is  important  to  note  that  this  selection  is  not  without  compromise.  The  most 
significant  limitation  of  the  trial  functions  is  with  regard  to  the  horizontal  deformation 
u.  Physically,  it  is  expected  the  w(.rj,')  behave  as  an  odd  function  as  discussed  above. 
However,  in  this  selection  of  trial  function,  a  positive  value  of  w(.vj-)  exists  at  the  origin. 
This  is  contrary  to  the  physical  behavior  of  the  problem  and  will  lead  to  some  error. 
However,  considering  that  the  magnitude  of  this  deformation  in  the  elastic  range  is 
small,  this  error  is  believed  to  be  limited.  Another  consequence  of  this  compromise  is 
the  existence  of  non-zero  shear  strain  at  jc  =  0. 

Another  less  severe  limitation  is  a  restriction  on  the  order  of  the  exponents  in 
order  to  maintain  zero  stress  at  the  boundaries.  Since  a,  and  ay  are  combinations  of  ex 
and  iy,  both  normal  strains  must  be  zero  at  the  boundaries  to  ensure  that  stresses  are 
zero  at  these  locations.  Since  the  normal  strains  are  first  derivatives,  this  requires  that 
b  and  c  are  at  least  equal  to  2. 

It  is  worthwhile  to  note  that  most  of  the  considerations  discussed  above  far  ex- 
ceed the  requirements  stipulated  by  the  Rayleigh-Ritz  method.  The  objective  has  been 
to  utilize  trial  functions  that  closely  match  the  physical  nature  of  the  problem  in  an  ef- 
fort to  maximize  accuracy. 

The  specific  trial  functions  used  in  this  study  are  listed  below.  The  horizontal 
deformation  was  assumed  to  be; 
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rri 

u(Xxy)  =  Yan(H-yf+n\  \  -  x)2  (3.3; 


n=\ 


where  //  and  L  represent  the  height  and  length  of  the  bearing  foundation,  respectively. 
Summation  was  done  for  n  equal  1  and  4.   The  vertical  deformation  was  assumed  to  be; 

m 

v{xy)  =  Ybn(H-yf(  -j  -  xf+n).  (3.4) 

1=1 

As  discussed  previously,  manipulation  of  the  exponents  enables  the  trial  function  results 
to  be  matched  with  the  analytical  solution.  As  will  be  illustrated  in  the  results,  the  ex- 
ponents chosen  in  the  above  functions  achieve  this  goal  sufficiently. 

With  the  deformations  chosen,  the  stresses  and  strains  can  be  calculated  for  use 
in  the  calculation  of  the  foundation  strain  enersv.    These  values  are  shown  below: 


tx  =  Y-2an(H-y){*+r,)(-^-x)  (3. 


n=\ 


n 

Zy  =  Yj-2bn(H-y)(j--xf+n)  (3.6) 


1(\-v)ex+vev]  (3.7a) 


(1  +v)(l-2v) 


a^(l+v,fl-2v)C(1-V)^+V^] 


(3.7b) 


where  E  and  v  are  Young's  modulus  and  Poisson's  ratio,  respectively.     For  the  shear 
stress  and  strain, 


m 


|t.y-(4  +  .K(II-/«(frf  (3.S) 


oy       i—i 


JL  =  Yj-(l+n)bn(H-y)2(j--xf  (3.9) 
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y       -f.  +  f.  (3.10) 

y       cy        ox 

Txy=GvXy  (3.11) 

where  G  is  the  shear  modulus. 

Using  the  above  quantities,  the  strain  energy  U  can  be  calculated.  From  the  definition 
of  strain  energy  applied  to  a  two  dimension  deformation  field, 


■\\ 


ii 
—  ( axix  +  oyzy  +  zxyyxy)dydx  (3.12) 


o 


Because  of  symmetry  about  the  origin,  strain  energy  can  be  calculated  for  half  of  the 
domain  and  doubled.  With  strain  energy  calculated,  the  total  potential  for  the  system 
can  be  found  from. 

n  =  U-FS  (3.13) 

where, 

F=  external  force  per  unit  length  applied  to  the  roller 
(3  =  vertical  distance  traveled  by  the  roller. 

The  quantity  FS  represents  the  work  done  by  the  roller  on  the  bearing  foundation. 

B.     FINITE  ELEMENT  APPROACH 
1.     Total  Potential  Derivation 

The  finite  element  method  can  be  employed  to  solve  the  contact  problem.  A 
finite  element  mesh  can  be  constructed  to  approximate  the  behavior  of  an  elastic  foun- 
dation subjected  to  line  contact  loading  from  an  rigid  roller.  The  resultant  interaction 
between  the  foundation  and  roller  enables  the  calculation  of  the  foundation's  strain  en- 
ergy and  the  subsequent  calculation  of  the  total  potential  energy.  By  again  utilizing  the 
optimization  techniques  discussed  in  Chapter  II,  the  equilibrium  position  can  be  deter- 
mined.  Thus  the  contact  stresses  can  be  calculated  throughout  the  body. 

The  objective  is  to  derive  a  means  of  calculating  the  total  potential  of  the  system 
by  application  of  the  finite  element  technique.    Total  potential  energy  is  defined  as, 
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n  =  U  -  F5 


(3.14) 


where, 

U=  internal  strain  energy 

F=  external  force  per  unit  length  applied  to  the  roller 

J  =  vertical  distance  traveled  by  the  roller. 

The  strain  energy  of  the  system  can  be  found  from. 


U  = 


—  \oxzx  +  oyey  +  rxyyxy)JQ 


3.15) 


where  Q  represents  the  problem  domain.    This  can  be  expressed  in  matrix  form  as, 


-T{c}T{a}ciQ 


Q 


3.16) 


where. 


l£)     -  icx   cy    "<'xy} 


{a}    =  [ox  a    t    } 


On  the  element  level. 


U  = 


±{e}T{a}dQ 


Q! 


(3.17) 


where  i  represents  analysis  of  the  i'h  element. 


The  stress  matrix  can  be  expressed  as, 


{*}  =  [£]{£} 


(3.18) 


where  [D]  represents  the  material  property  matnx. 
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For  a  condition  of  plane  strain,  the  stiffness  matrix  can  be  expressed  as, 


[0] 


For  a  plane  stress  condition. 


E 


1  +v)(l 


2v] 


1   —  V  V 

V  1  —  V 

0  0 


0 

0 

1  -  2v 


(3.19a) 


[D]  = 


(1-v' 


1 

V 

0 

V 

1 

0 

0 

0 

1  -  V 

3.1%) 


The  development  of  this  technique  will  use  linear  triangular  elements.  The  method, 
however,  can  be  applied  to  any  type  element.  For  linear  triangular  element,  the  defor- 
mations take  the  followin2  form: 


u  =  Hiii\  +  H-^ih  +  H-u 


i"i 


3"3 


v  =  //.v,  +  H,\2  +  H,\ 


,»,  -r  111*2  ^  M3'3 


where  the  shape  functions  H,  are  defined  as: 


(3.20) 
(3.21) 


//, 


=  C(-v^'3  ~  x3>'2)  +  CV'2  ~}'3)x  +  (x3  ~  *:)>']  ~ 


(3. 


■>->} 


H2  =  [(-ViXi  -  xy3)  +  CV'3  -}'\)x  +  (x\  ~  x3]y]  - 


IA 


(3-23) 


Hi  =  [(*i>2  ~xiy\)  +  Oi  -y-d*  +  (*2  -  *ilv]  37 


where, 

x(,y,  =  coordinates  for  node  i 
A  =  element  area   [Ref.  14  ]. 

The  strain  matrix  can  be  expressed  as, 


(3.24) 
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w.= 


4-  o 

ex 

o    f 

g        c 
c  v      dx 


Substituting  equations  (3.20)  and  (3.21)  into  equation  (3.25)  yields, 


o//, 

0 

6H2 

0 

c//3 

0 

dx 

ex 

0 

cH{ 

0 

a//2 

0 

c//3 

°y 

cy 

cy 

c//, 

cHx 

c//: 

o//2 

BH3 

dH3 

c    = 


cy        ex        cy        ex         cy        ex 
In  abbreviated  notation,  equation  (3.26)  is  expressed  as, 

For  the  linear  triangular  element,  [5]  reduces  to, 


r         -i 

"l 

Vl 

u2 

v2 

''3 

V3 

[Z?]  =  17 


>2->'3         0         >'3-yi  0         y,  -y2         0 

0  ,v3  -  x2  0  .V!  -  ,v3  0  x2  -  -V, 

-r3  -  x2     }'2  -}'2      x\  -  x2     }'3  ~}'\      x2  ~  x\     }'\  ~}'2 


i.28) 


Returning  to  the  element  strain  energy  calculation,  equation  (3.17), 


U  = 


& 


±{c}T{o}d£l 


and  substituting  (3.18)  into  (3.17), 


i 


±-{,}\D^z}dV. 


(3.29) 
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Substituting  equation  (3.27)  yields, 


U  = 


\{lB-]{d])\D^B^{d)dV 


q: 


r  =  4-{^r[5]r[D][5]{c/}f  d\ 


l-  =  -T{cf}T[Bf[D]iBl4i{J}  (3.30) 


where  i  =  unit  depth.    Defining  the  element  stillness  matrix  [A.-] 


[A1  =  [*]r[D][5].-!/  13.31) 

then  the  strain  energy  per  element  equals, 

U  =  ±-{d}tK\{d\  (3.32) 

The  element  stiiTness  matrix  can  be  expanded  into  the  global  stillness  matrix. 
With  strain  energy  now  determined,  total  potential  can  be  determined  from  equation 
(3.14). 

2.     Optimization  and  Static  Condensation 

As  with  the  previously  developed  model,  the  augmented  Lagrange  multiplier 
method  will  be  utilized  to  determine  the  equilibrium  condition  via  the  theorem  of  mini- 
mum potential.  The  objective  function  is  again  the  total  potential.  In  this  case,  how- 
ever, the  design  variables  are  the  nodal  deformations,  ut  and  v„  for  non-fixed  nodes. 
Constraint  equations  are  developed  in  a  similar  manner  as  discussed  in  Chapter  II  to: 
ensure  that  one  body  does  not  violate  space  occupied  by  the  other. 

Since  the  nodal  deformations  are  represented  as  the  optimization  design  vari- 
ables, the  number  of  design  variables  will  equal  twice  the  number  of  non-fixed  nodes. 
For  a  simple  mesh,  a  direct  application  of  this  procedure  will  likely  yield  accurate  results. 
However,  it  is  known  that  the  accuracy  of  the  optimization  routine  declines  as  the 
number  of  design  variables  increases.  Hence,  for  complicated  finite  element  meshes, 
solution  accuracy  will  be  adversely  affected  by  the  large  number  of  design  variables. 
Therefore  a  procedure  must  be  adopted  to  eliminate  the  need  for  assigning  design  vari- 
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ables  to  nodes  where  information  is  not  necessary  for  evaluating  a  solution.    This  pro- 
cedure is  known  as  static  condensation. 

Static  condensation  has  been  utilized  by  References  2  and  3  in  an  effort  to  re- 
duce computer  computational  time.  The  idea  behind  static  condensation  is  the  reor- 
ganization of  the  global  stiffness  matrix.   A  finite  element  problem  can  be  expressed  as, 


[KJ{u}  =  {F} 


(3.33) 


where. 

[K,]  =  the  stiffness  matrix 
{«}  =  the  deformation  vector 
{F}  =  the  force  vector 

It  is  desired  to  reorganize  this  system  of  equations  into  the  following: 


^11     ^12 
Kt,     Ao-> 


n 

of 


(3.34) 


The  vector  m,  contains  the  essential  nodes  while  vector  u2  contains  non-essential  nodes. 
Essential  nodes  are  those  where  boundary  conditions  are  applied  and  nodes  that  are  as- 
signed optimization  design  variables.    By  matrix  manipulation. 


{u2)  =  -lK22TXK2l-]{ux}. 


(3.35) 


Therefore,  the  displacement  vector  can  be  expressed  as, 


u2}"\-LK22TlLK2l2M] 


[/] 


-  IK22T\K2{\ 


(3.36) 


where  I  represents  the  identity  matrix. 


Substituting  this  equation  into  the  global  counterpart  of  equation  (3.32), 


,;  1       f         ,71 

V  =y  (Mi) 


IK22Y\K2{] 


\KS~] 


[^22]     [^21] 


fo 


(3.37) 


By  defining  the  reduced  stiffness  matrix  [A-],  the  above  reduces  to. 
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U  =  ±-{u{}T[K}{Ul)  (3.38) 

As  discussed  at  the  beginning  of  this  section,  the  ADS  design  variables  are  as- 
signed as  the  horizontal  and  vertical  deformations  at  all  non-fixed  nodes.    With  the  in- 
tegration of  static  condensation,  design  variable  assignments  are  further  restricted  ten 
non-fixed,  non-condensed  nodes.   The  procedure  is  now  in  place  for  calculation  of  strain: 
energy  and  total  potential  energy.  In  summary,  a  flow  chart  of  the  solution  procedure' 
utilized  in  this  chapter  is  shown  in  figure  9. 
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(     Initialize  Program 

V I y 


Define  Optimization  Design  Variables 


Calculate  Element  Stiffness  Matrix  [K] 


Assemble  Global  Stiffness  Matrix  (KG] 


Reduce  Degree  of  Freedom  by  Static  Condensation  [ KG ] 


Alter  Design 
Variables 


Assign  Design  Variables  to  (d) 


Calculate  Strain  Energy  U 
U={d}t(KG]{d) 


Calculate  Total  Potential 


Define  Constraints 


Call  Optimizer  to  Minimize 
Total  Potential 


Post-Process 


C   End  ) 


Figure  9.      Finite  element  method  applied  to  contact  problem 
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IV.     RESULTS  AND  DISCUSSION 

A.     PROCEDURE  VALIDATIONS 
1.     Rayleigh-Ritz  Method  Results 

In  Chapter  III  an  approximation  technique  was  developed  to  solve  the  contact! 
problem  via  the  Rayleigh-Ritz  method.    As  discussed,  two  trial  functions  that  approxi- 
mated the  deformation  field  were  selected  in  terms  of  unknown  coefficients.    The  hori- 
zontal deformation  was  assumed  to  be; 

m 
n=] 

The  vertical  deformation  was  assumed  to  be; 

m 

vixvv)  =  Yjbn(H-y)2(-T-x)0+n) 

n=\ 

Using  the  above  deformation  fields,  theory  of  elasticity  relationships  and  the  definition 
of  strain  energy  were  employed  to  obtain  an  expression  for  the  total  potential  energy  of 
the  system  shown  in  Figure  1.  Numerical  minimization  techniques  were  then  employed 
to  determine  the  equilibrium  condition  and  the  contact  stresses. 

To  illustrate  the  application  of  this  method,  an  isotropic  material  with  the  fol- 
lowing properties  was  selected: 


E  =  200  GPa 

v  =  0.3 

G  =  76.9  GPa 


As  stated  in  Chapter  II,  this  problem  was  selected  for  development  of  this  technique 
because  an  analytical  solution  is  available  as  a  result  of  the  work  done  be  Hertz.  It  is 
desired  to  use  this  analytical  solution  to  choose  a  roller  size,  load,  and  problem  domain 
that  can  be  used  to  accurately  simulates  the  contact  phenomenon. 

Contact  stresses  as  well  as  the  size  of  the  surrounding  region  of  influence  are 
strongly  affected  by  the  size  of  the  contact  zone.    Naturally,  as  the  size  of  the  contact 
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zone  increases,  the  load  is  distributed  over  a  larger  area  and  contact  stresses  decrease. 
The  extent  of  the  affected  subsurface  zone  also  decreases.  From  the  analytical  solution, 
the  contact  zone  size  is  defined  by  the  externally  applied  force,  the  material  properties 
and  the  diameter  of  the  roller  [  Ref.  13  ].  Hence  for  a  given  material,  the  roller  size  and 
the  external  load  define  the  contact  zone  size  and  the  resulting  stresses. 

Figure  2  defined  the  width  of  the  contact  zone  as  2a.  Figure  3  shows  the  ana- 
lytical solution  of  the  contact  problem.  This  figure  shows  the  decrease  of  the  normal 
stresses  as  a  function  of  half-contact  zones  (a)  awa\  from  the  contact  point.  As  shown. 
ay  decreases  more  gradually  than  o,  .  Therefore  the  decay  o^  ay  is  the  limiting  factor  in 
defining  the  domain  beyond  which  strain  energy  contributions  are  negligible.  From 
Figure  3.  it  is  estimated  that  the  contact  phenomenon  can  be  accurately  modeled  by 
examining  a  region  equal  to  approximately  five  half-zones  (5a). 

Since  a  numerical  integration  technique  was  used  to  perform  the  double  inte- 
gration required  by  Equation  3.12.  the  dimensions  were  selected  for  numerical  conven- 
ience. Referring  to  Figure  8.  height  H  and  length  L  were  selected  as  1  and  2  meters, 
respectively.  Due  to  the  problem's  symmetry,  half  the  foundation  was  analyzed.  This 
enabled  the  double  integration  to  be  conducted  between  the  limits  of  0  and  1.  Since  the 
foundation  height  H  has  been  set  to  1  m,  it  is  desired  to  have  this  distance  equal  to  5 
contact  zones  (5a)  as  described  above. 

Using  the  analytical  solution,  a  load  and  roller  diameter  were  selected  that  cre- 
ated a  contact  zone  such  that  the  foundation  height  H  was  equal  to  a  distance  of  5a. 
The  only  additional  restriction  was  that  the  resulting  contact  stresses  remained  within 
the  elastic  range  of  the  material.  Yield  stress  was  assumed  to  be  300  MPa.  The  load 
and  roller  radius  combination  used  in  this  study  are; 

Load:     90  MN 

Radius:  75  m 

The  values  were  obtained  using  the  analytical  solutions  found  in  Reference  13.  The 
latter  dimension  is  unrealistic.  However,  as  stated  above,  it  is  a  result  of  the  selection 
of  base  dimensions  in  the  interest  of  numerical  convenience.  This  model  simply  repres- 
ents a  scaled-up  examination  of  the  base  material  in  a  small  region  adjacent  to  the  con- 
tact zone. 

Comparisons  of  the  approximated  contact  stresses  with  the  analytical  solution 
along  the  axis  of  symmetry  are  shown  in  Figures  10  and  11.   Both  figures  are  normalized 
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graphs  of  the  stress  distribution.  As  shown  in  Figure  3,  the  maximum  stress  occurs  at 
the  point  of  contact.  As  shown,  this  stress  is  equal  to  ax  and  oy  at  the  point  of  contact: 
Figure  10  is  a  comparison  for  ay  .  This  figure  shows  a  stress  distribution  that  closely 
approximates  the  analytical  solution.  Figure  1 1  shows  the  analytical  and  approximate 
stress  distribution  for  ax.  It  is  apparent  that  this  method  over  approximates  this  stress.- 
As  stated  in  Chapter  III,  one  of  the  benefits  of  this  trial  function  is  the  ability  to  change 
the  exponents  to  match  analytical  solutions  with  approximate  solutions.  A  brief  expla- 
nation of  the  choice  of  exponents  used  in  this  solution  technique  follows. 

While  selecting  exponents,  it  must  be  understood  that  ax  and  ay  are  both  func- 
tions of  horizontal  and  vertical  deformations.  Therefore  changing  the  exponent  of  one 
deformation  to  alter  the  stress  in  one  direction  will  influence  the  behavior  of  the  other. 
In  the  case  of  Figure  11,  it  appears  that  a  modification  is  required.  An  increase  in  the 
exponent  of  the  y-portion  of  the  horizontal  deformation  function  seems  appropriate  to 
increase  the  rate  at  which  ax  decays.  However,  this  action  will  have  the  undesirable  ef- 
fect of  decreasing  the  rate  of  decay  of  ay.  It  is  possible  to  counter  by  decreasing  the 
exponents  of  the  vertical  deformation.  However,  as  stated  in  Chapter  III.  in  order  to 
maintain  a  zero  stress  boundary  condition  at  x=L  2  and  y=H  the  exponents  of  all 
terms  must  be  greater  than  or  equal  to  2.  Thus  a  compromise  must  be  reached.  It  is 
believed  that  since  ay  is  the  dominant  term,  priority  should  be  placed  on  ensuring  ay  is 
as  close  as  possible  to  the  analytical  solution.  Accordingly,  a  decision  was  made  to  ac- 
cept the  over  estimation  of  ax  as  shown  in  Figure  11.  In  this  case,  the  estimation  of  ax 
is  conservative. 

The  contours  of  the  Rayleigh-Ritz  solution  are  shown  in  Figures  12  through  15. 
It  has  been  determined  that  this  method  can  accurately  predict  contact  stresses  resulting! 
from  line  contact  between  a  roller  and  flat  plane.  Table  1  compares  the  approximated! 
contact  stresses  and  the  analytical  results  at  the  point  of  contact. 


Table   1.     RAYLEIGH-RITZ  RESULTS  AT  CONTACT  POINT 


Current  Model 

Analytical  Solution 

a,  (MPa) 

285.6 

289.7 

ov  (MPa) 

301.9 

289. 7 
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Figure   12.      Stress  contour  of  «ry  from  Rayleigh-Ritz  method 

2.     Finite  Element  Method  Results 

In  Chapter  III,  an  approximation  method  was  developed  to  solve  the  contact 
problem  using  the  finite  element  method.  A  method  of  calculating  a  system's  strain  en- 
ergy and  the  total  potential  energy  was  investigated.  In  addition,  the  use  of  static 
condensation  to  improve  optimization  efficiency  was  described.  As  discussed,  the  nu- 
merical minimization  techniques  were  again  utilized  to  determine  the  equilibrium  posi- 
tion. A  means  of  employing  these  techniques  to  evaluate  the  contact  phenomenon  was 
introduced.  In  this  section,  a  simple  contact  problem  will  be  investigated  to  validate  the 
algorithms  used  to  calculate  strain  energy  and  those  used  to  implement  static 
condensation.  With  confidence  in  these  algorithms,  the  more  complex  roller-foundation 
problem  will  be  examined. 

a.      Two  Thin  Plates  in  Contact 

The  procedure  developed  was  first  validated  on  a  simple  contact  problem 
the  solution  of  which  was  known.    In  this  example,  two  thin  bodies  in  plane  stress  were 
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Figure   13.      Stress  contour  of  <x,  from  Rayleigh-Ritz  method 

studied.  As  shown  in  Figure  16,  one  body,  restrained  on  one  edge  and  subjected  to  a 
horizontal  load,  comes  in  contact  with  a  second  body  rigidly  supported  on  three  sides. 
The  Finite  element  model  developed  to  solve  this  problem  is  composed  of  14  linear  tri- 
angular elements  as  shown  in  Figure  17.  The  objective  of  the  fortran  program  developed 
to  solve  this  problem  was  to  calculate  the  total  potential  energy  of  the  system  using  the 
finite  element  technique  and  the  method  of  static  condensation.  With  this  accomplished, 
the  equilibrium  position  can  then  be  found  via  the  augmented  Lagrange  multiplier 
method. 

The  objective  function  for  this  problem  is  the  total  potential  energy, 
Equation  3.14.  Referring  to  Figure  16  and  17  the  constraints  imposed  on  the  system  are 
expressed  as: 

u&  <  0.005  +  uu 
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Figure   14.      Strain  contour  of  ty  from  Rayleigh-Ritz  method 

tig  <  0.005  +  ul2 
where  u  represents  the  horizontal  deformation  of  node  /. 

As  shown  in  Figure  17,  seventeen  nodes  were  used  to  model  the  system. 
Each  node  has  been  assigned  horizontal  and  vertical  deformation  variables.  Accord- 
ingly, the  degree  of  freedom  for  this  system  is  twice  the  number  of  nodes.  As  discussed 
in  Chapter  III,  static  condensation  requires  the  identification  of  essential  nodes  and 
non-essential  nodes.  Essential  nodes  are  those  nodes  where  ADS  design  variables  are 
assigned  and  boundary  conditions  are  applied.  Referring  to  Figure  17,  node  3  is  the 
point  of  load  application  and  must  be  assigned  a  design  variable.  Nodes  8,  9,  11,  and 
12  are  assigned  design  variables  in  order  to  define  the  constraint  equations  described 
above.  After  eliminating  all  fixed  nodes  from  consideration,  nodes  2,  5,  and  6  are  the 
only  candidates  for  condensation.    As  discussed  in  Chapter  III,  the  global  stillness  ma- 
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Figure   15.      Strain  Contour  et  from  Rayleigh-Ritz  method 

trix  is  rearranged  according  to  Equation  3.34.  For  this  problem,  the  vector  containing 
the  non-essential  nodal  information.  u2,  is  arranged  as  follows; 

{w2}r={w2    v2    u5    v5    u6    v6} 

where  u,  and  v\  represent  the  horizontal  and  vertical  deformation  of  node  /,  respectively. 
Strain  energy  and  total  potential  energy  were  calculated  according  to  Equations  3.38  and 
3.14.  The  latter  was  minimized  using  the  augmented  Lagrange  multiplier  method  de- 
scribed in  Chapter  II. 

The  solution  obtained  from  this  simple  problem  were  compared  with  the 
results  obtained  from  Y.W.  Kvvon  and  J.E.  Akin  [  Ref.  15  ]  and  are  shown  in  Table  2. 
The  solutions  were  in  agreement.    It  was  concluded  that  a  satisfactory  procedure  was  in 
place  to  solve  the  more  complex  roller-foundation  problem. 
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Figure   16.      Tno  thin  plates  in  contact 


Table  2.     FINITE  ELEMENT  RESULTS:  TWO  PLATES  IN  CONTACT 


Load  (N) 

Deformation 

Current  Model 

Reference  15 

1.0x10s 

No  Contact 

No  Contact 

1.0x10° 

"9 

"12 

.503xlO-2 
.336X10"4 

.503x10-' 

.319x10^ 

l.OxlO1 

"9 

un 

.702xlO-2 
.201xlO~2 

.734x10^ 
.235xlO"2 

b.     Roller- Foundation  Contact  Problem 

A  finite  element  grid  composed  of  512  linear  triangular  elements  was  con- 
structed to  model  the  roller-foundation  assembly  shown  in  Figure  1.  Because  of  the 
symmetry  of  the  problem,  one-half  of  the  foundation  was  modeled.  A  refined  mesh  was 
constructed  in  the  vicinity  of  the  point  of  contact.  The  mesh  is  shown  in  Figure  18 
where  the  origin  represents  the  point  of  contact.    The  domain  dimensions  are  similar  to 
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Figure   17.      Finite  element  mesh  for  two  plates  in  contact 

those  chosen  in  the  Rayleigh-Ritz  method  discussed  in  Part  A  of  this  chapter.  The  roller 
radius  was  chosen  as  75  meters  and  the  half-domain  dimensions  are  1.40  x  1.40  meters,  j 
As  discussed  in  Part  A,  these  dimensions  represent  an  analysis  of  the  region  immediately 
adjacent  to  the  contact  zone  and  are  a  result  of  the  local  nature  of  the  contact  problem. 
Constraint  equations  were  constructed  according  to  the  discussion  of 
Chapter  II  Part  B.    Boundary'  conditions  were  imposed  in  the  following  manner: 

•  Horizontal  and  vertical  deformations  were  prohibited  on  the  remote  mesh  bound- 
aries (i.e.,  m(  1.40, y)  =  v(  1.40, y)  =  0  and  u(jr,1.40)  =  v<x,1.40)  =  0). 

•  Horizontal  deformation  was  prohibited  on  the  axis  of  symmetry  (i.e.,  u(0,y)  =  0). 
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Figure   18.       Finite  element  mesh  for  roller-foundation  problem 
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Since  the  system  has  289  nodes,  the  resulting  578  degrees  of  freedom  ne-1 
cessitated  the  utilization  of  static  condensation.     Referring  to  the  contact  surface,  the 
deformation  variables  that  correspond  to  nodes  on  this  surface  are  required  for  con- 
straint equations.    The  nodes  that  comprise  the  other  three  borders  of  the  domain  arc: 
all  subject  to  boundary  conditions.   Consequently,  the  interior  nodes  of  the  mesh  are  the 
nodes  that  are  candidates  for  static  condensation.    In  this  model,  all  interior  nodes  were 
condensed.   The  original  system  was  reduced  from  578  to  128  degrees  of  freedom.   After; 
application  of  the  boundary  conditions,  there  were  46  possible  deformations,  a  suffi- 
ciently small  number  of  design  variables  for  the  optimization  algorithm.   One  additional 
design  variable  was  used  to  represent  the  distance  of  travel  by  the  roller.    This  value  is' 
needed  to  calculate  the  work  done  by  the  roller  on  the  bearing  foundation. 

As  before.  Equation  3.38  was  used  to  calculate  the  system's  strain  energy. 
Following  calculation  and  the  subsequent  minimization  of  the  total  potential  energy,  a 
post-processing  procedure  was  followed  to  determine  the  contact  stresses.  The  output 
of  the  optimization  routine  represents  the  nonzero  components  of  the  {«,}  vector.  In 
order  to  calculate  stresses  throughout  the  body,  the  remaining  deformations  contained 
within  the  condensed  vector  {;/,}  must  be  determined.  This  vector  is  calculated  directly 
using  Equation  3.35.  With  deformations  known  throughout  the  domain,  strains  can  be 
determined  by  applying  Equation  3.27  to  each  element.  The  subsequent  application  of 
Equation  3.18  enables  determination  of  stress  for  each  element. 

To  illustrate  the  capability  of  this  method,  an  isotropic  material  with  the 
following  properties  was  chosen: 

E  =  240  GPa 
v  =  0.3 
G  =  92.3  GPa 
Load  =  90.0  MPa 

Figure  19  shows  the  deformation  resulting  from  the  loading.  For  clarity,  the  defor- 
mations have  been  magnified  100  times  their  original  values.  Comparisons  of  the  stress 
distributions  with  the  analytical  solution  along  the  axis  of  symmetry  are  shown  in  Fig- 
ures 20  and  21.  These  figures  are  similar  to  Figures  10  and  1 1  and  represent  normalized 
versions  of  contact  stresses.  As  shown  in  these  figures,  this  method  is  a  good  approxi- 
mation of  the  stress  distribution  in  the  foundation  of  a  loaded  roller  bearing.   If  the  mesh 
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was  more  refined  near  the  contact  zone  and  the  domain  extended  further,  the  agreement 
between  the  numerical  and  analytical  solutions  would  be  better,  figures  22  and  23  rep- 
resent normal  stress  contours  of  this  problem.  Figures  24  and  25  show  normal  strains. 
Table  3  shows  a  comparison  of  the  results  of  this  model  and  the  analytical  solution  at 
a  selected  element  in  the  region  of  contact. 

Table  3.     COMPARISON  OF  STRESSES  NEAR  THE  POINT  OF  CONTACT 


FEM  Solution 

Analytical  Sol- 
ution 

a,  (MPa)  at  x  =  0.0137.  y  =  0.0273 

281.3 

272.0 

ex,  (MPa)  at  x  =  0.0137,  y  =  0.0273 

315.7 

313.8 

B.     APPLICATION 

The  preceding  section  illustrated  that  the  contact  problem  can  be  accurately  simu- 
lated using  the  methods  developed  in  Chapter  III.  It  is  the  objective  of  this  section  to 
show  how  this  approach  can  be  applied  to  a  contact  problem  in  a  composite  plate  sub- 
jected to  low-velocity  impact. 

A  multi-ply  laminate  model  has  been  constructed  to  investigate  the  response  of 
composite  materials  to  low  velocity  impact.  It  has  been  found  that  composite  bodies 
subject  to  impact  damage  commonly  fail  due  to  delamination.  Sandwich  composites  are 
currently  being  considered  for  use  as  turbine  blades.  It  would  be  beneficial  to  acquire 
an  understanding  of  the  behavior  of  sandwich  materials  to  impact  damage. 

In  order  to  accomplish  this  task,  a  clamped  composite  beam  similar  to  the  one  de- 
picted in  Figure  26  has  been  modeled.  The  beam  length  is  25  cm.  Beam  thickness  is  2.5 
cm.  Because  of  symmetry,  half  the  beam  was  modeled  with  256  bilinear  elements.  As 
shown  in  Figure  27,  the  mesh  is  refined  near  the  point  of  contact,  the  origin  of  the  mesh. 

The  major  assumptions  of  this  model  are  that  the  beam  is  in  a  condition  of  plane 
strain  and  that  the  dynamic  effect  of  the  impact  can  be  neglected.  Reference  10  ap- 
proximated the  loading  resulting  from  low  velocity  impact  as  a  sinusoid.  In  this  study, 
the  peak  load  will  be  examined  to  study  the  maximum  normal  and  shear  stresses.  Ac- 
cordingly, the  stress  distributions  obtained  from  this  study  will  represent  a  'snap-shot' 
in  time  of  the  response  of  the  body  at  the  instant  of  maximum  loading. 

The  finite  element  program  developed  to  solve  this  problem  is  sufficiently  flexible  to 
alter  the  material  stiffness  matrix  [D]  shown  in  Equation  3.18  during  construction  of  the 
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Figure   19.      Deformed  finite  element  solution:  deformation  magnified  100  times 
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Figure  20.      Analytical  solution  vs.  approximate  ay  from  finite  element  results 
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Figure  21.      Analytical  solution  vs.  approximate  <xx  from  finite  element  results 
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Figure  22.      Stress  contour  a   from  finite  element  results  (increment  0.01  GPa) 
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Figure  23.      Stress  contour  <x,  from  finite  element  results  (increment  0.02  GPa) 
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Figure  24. 


Strain  contour  sy  from  finite  element  results  (increment  0.000035) 
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Figure  25.      Strain  contour  £,  from  finite  element  results  (increment  0.00001) 
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Figure  26.      Clamped  composite  beam 
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Figure  27.      Finite  element  mesh  for  clamped  beam  model 

finite  element  global  stiffness  matrix.  Therefore  by  defining  the  layup  for  the  laminate, 
the  lamina  stiffness  matrices  can  be  varied  from  element  to  element  to  accurately  model 
the  behavior  of  the  body.  This  enables  a  variety  of  laminate  layups  and  lamina  thick- 
nesses to  be  studied. 

The  sandwich  materials  used  in  this  study  are  composed  of  an  isotropic  interior 
material  and  orthotropic  exterior  laminae.  Since  a  condition  of  plane  strain  was  as- 
sumed, the  isotropic  material  stiffness  matrix  is  given  by  Equation  3.19a.  For  the 
orthotropic  exterior  laminae,  the  material  stiffness  matrix  is  given  by  Equation  4.1 
[Ref.  16  ]. 
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where. 
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v::v:,  -  V3VJ,-  v31v13-  2l\,v3:v13 


r= 

£,£2£3 

£  =  Young's  Modulus  in  /"  direction 

v .  =  Poisson's  ratio  for  lateral  contraction  in/'"  direction  resulting  from  loading  in 
the  r*  direction 

The  derivation  of  the  total  potential  energy  calculation  in  Chapter  III.  Part  B  was 
done  using  linear  triangular  elements.  Since  bilinear  elements  were  used  in  this  model, 
calculation  of  the  element  stiffness  matrix  was  more  computationally  intensive.  Indi- 
vidual entries  of  the  element  stiffness  matrix  were  obtained  from  Reference  17.  Other- 
wise, the  calculation  of  total  potential  energy  was  identical  to  the  procedure  outlined  in 
Chapter  III. 

There  were  two  groups  of  boundary  conditions  applied  to  the  problem.  Along  the 
clamped  edge,  deformation  was  prohibited.  In  addition,  horizontal  deformation  was 
prohibited  along  the  axis  of  symmetry.  As  before,  the  application  of  static  condensation 
requires  the  identification  of  essential  nodes  the  information  of  which  is  contained  within 
the  {w,}  vector.  In  addition  to  the  essential  nodes  associated  with  the  above  boundary- 
conditions,  the  nodes  along  the  contact  surface  are  needed  for  the  constraint  equations. 
All  other  nodes  were  condensed  out. 

To  illustrate  this  application  of  problem  solving,  sandwich  material  with  the  fol- 
lowing properties  were  studied: 

Exterior  Laminae 

£u  =   170  GPa 

Ex  =   11.8  GPa 

Gu  =  5.2  GPa 

vi2  =  0.33 
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Isotropic  Core 
E  =  2.24  GPa 
G  =  0.84  GPa 
v  =  0.35 

The  beam  was  loaded  by  contact  with  a  10  cm  radius  ball.  The  transmitted  force 
was  250  N.  Some  unexpected  trends  were  observed  in  the  equilibrium  position  deter-| 
mined  by  the  optimizer.  By  examining  the  deformations  along  the  axis  of  symmetry,  a 
gradually  decreasing  trend  in  deformation  moving  away  from  the  point  of  contact  was 
interrupted  in  lamina  of  significantly  decreased  stiffness.  It  is  believed  that  this  difficulty 
resulted  from  an  inability  of  the  optimization  routine  to  approximate  deformations 
through  regions  containing  very  different  orders  of  strain  energy.  In  spite  of  these  dif- 
ficulties, some  critical  information  was  obtained  from  this  program.  As  discussed  in  the 
beginning  of  this  study,  one  of  the  greatest  difficulties  of  the  contact  problem  is  the  de- 
termination of  the  size  of  the  contact  zone.  Fortunately,  the  size  of  the  contact  zone  can 
be  readily  determined  by  examining  the  output  from  the  constraint  equations.  By  com- 
paring the  ball  radius  (r)  and  the  distance  between  the  ball  center  and  the  node  (r-),  it 
can  be  determined  if  a  node  is  in  contact.  This  condition  is  illustrated  in  Figure  28.  The 
distance  r'  to  the  im  node  is  given  by  the  equation: 


r'  =  v  (r  -  <5  +  vf  +  x? 


If  r'  is  greater  than  r,  the  node  is  not  in  contact  with  the  body. 

With  the  extent  of  the  contact  zone  known,  the  solution  to  this  problem  was  ob- 
tained by  applying  the  contact  boundary  condition  to  a  direct  finite  element  program. 
Since  the  validity  of  the  optimization  program  was  in  question,  this  method  was  applied 
using  a  0-90-0  layup  similar  to  one  used  in  a  study  conducted  by  Sun  and  Rechak  [Ref. 
10].  The  solution  obtained  from  the  current  approach  was  very  close  to  the  other  re- 
sults. 

With  confidence  in  the  procedure,  it  was  desired  to  examine  the  behavior  of  this 
model  to  various  layups.  The  objective  was  to  illustrate  how  this  solution  technique  can 
be  used  for  meaningful  research.  The  study  conducted  by  Sun  and  Rechak  analyzed 
methods  of  reducing  the  likelihood  of  composite  failure  due  to  delamination.  Of  par- 
ticular concern  is  the  magnitude  of  shear  stress  distribution  and  tensile  stress  in  the  y- 
direction,  the  two  predominant  causes  of  delamination  failure.    Using  materials  with  the 
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Figure  30.      Stress  contour:  xxy  for  clamped  beam  model 

thickness  of  the  exterior  layer  increases,  the  maximum  value  of  this  stress  increases. 
However,  this  stress  is  always  compressive  at  the  interface  with  the  core.  By  comparing 
the  magnitudes  of  the  shear  and  transverse  normal  stresses  at  the  interface,  it  would 
appear  that  if  delamination  was  to  occur  at  this  interface,  it  is  more  likely  to  be  caused 
by  high  shear  stresses. 

Comparisons  of  the  bending  stresses  at  cross  section  A  are  shown  in  Figures  39,  40, 
and  41.  These  figures  show  that  as  the  thickness  of  the  exterior  layer  increases,  a  non- 
linear stress  distribution  intensifies  in  the  layer  closest  to  the  contact  surface.  This  trend 
would  indicate  that  beam  theory  is  unsuitable  for  estimating  bending  stress  through  this 
lamina.  This  nonlinear  behavior  is  local  to  the  contact  zone.  Figure  42  shows  the 
counterpart  for  Figure  41  at  cross  section  C.  The  stress  distribution  in  the  contact  layer 
is  approximately  linear.  Another  significant  observation  can  be  made  by  examining  the 
bending  stress  graphs.  As  the  thickness  of  the  exterior  layer  opposite  to  the  contact 
layer  increases,  the  stress  distribution  within  this  layer  transforms  from  purely  tensile 
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Figure  31.      Stress  contour:  <tx  for  clamped  beam  model 


behavior  to  compressivc-tensile  behavior.  This  is  significant  because  a  bending  crack 
initiated  by  tensile  stresses  tends  to  propagate  to  the  core  interface  and  cause  delami- 
nation.  The  presence  of  compressive  stress  within  this  layer  will  tend  to  slow  the  growth 
of  this  crack  toward  the  core. 

Thus  far  symmetric  layups  have  been  studied.  To  analyze  how  beams  with  asym- 
metric layups  respond  to  contact  loading,  two  beams  with  the  following  designations 
were  studied:  0(3)-ISO(6)-0(7)  and  0(8)-ISO(6)-0(2).  The  first  designated  layer  is  the 
lamina  closest  to  the  contact  surface.  Figures  43  and  44  show  the  shear  stress  distrib- 
utions for  these  two  layups.  As  before,  the  maximum  shear  stress  is  relatively  unaffected 
by  the  different  layups.  As  was  the  case  for  the  symmetric  beams,  a  thicker  exterior  layer 
close  to  the  contact  zone  results  in  a  more  gradual  transition  of  shear  stress  into  the 
core.  The  result  is  a  lower  shear  stress  at  the  interface  for  the  0(8)-ISO(6)-0(2)  case. 
As  seen  in  Figures  45  and  46,  the  trends  for  the  transverse  normal  stress,  ay,  were  the 
same    as    those    found    in    the    svmmetnc    beams.        Deflection    was    less    for    the 
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Figure  28.      Determination  of  nodal  contact 

properties  outlined  above,  a  sandwich  composite  beam  with  outer  fibers  aligned  to  0 
degrees  was  first  studied.  Since  the  finite  element  model  is  composed  of  16  layers,  the 
beams  studied  will  be  described  with  the  number  of  finite  element  layers  in  parenthesis 
following  the  layer  description.  For  example,  a  0(3)-ISO(  10)-0(3)  beam  is  composed  of 
10  isotropic  core  layers  within  3  layers  of  material  with  the  fibers  oriented  at  0  degrees 
on  the  top  and  bottom  of  the  beam. 

Three  symmetric  layups  of  varying  core  thickness  were  initially  studied.  The  beams 
have      the      following      designations:  0(3)-ISO(  10)-0(3),      0(4)-ISO(8)-0(4),      and 

0(5)-ISO(6)-0(5).  Figure  29  shows  the  deformed  0(3)-ISO(lO)-0(3)  beam  with  defor- 
mations magnified  100  times.  Stress  contours  for  this  beam  are  shown  in  Figures  30, 
31,  and  32.  Figure  30  shows  the  shear  stress  contour  for  the  loaded  condition.  This 
stress  is  of  particular  concern  since  delamination,  a  common  failure  mode  for  compos- 
ites, is  commonly  initiated  by  high  shear  stresses  or  tensile  transverse  normal  stresses. 
As  the  figure  shows,  a  very  high  stress  gradient  is  present  near  the  contact  zone.  As  the 
distance  along  the  beam  increases  away  from  the  contact  zone,  the  magnitude  of  the 
gradient  decreases  until  the  cross  sectional  shear  stress  distribution  becomes  parabolic. 
The  transverse  normal  stress  is  also  concentrated  around  the  contact  zone. 

Figures  33,  34,  and  35  show  the  cross  sectional  shear  stress  distributions  for  the 
three  symmetric  beams  described  above.    Three  separate  cross  sections  are  shown  on 


55 


s 


— r ' 7  S  10-° 

a'S  HWIZDKIfL  DISTBNDE  FROM  POINT  OF  CONTACT 


ia-5 


Figure  29.      Loaded  beam:  deformation  magnified  100  times 


each  oraph.  the  locat.ons  of  which  are  indicated  in  Ftgure  27.    The  "™*^^ 
onea'ch  gr  Ph  identify  the  lamma  interfaces.   These  figures  she,  that  a,,  stgntfi  an     c- 

<h«ir  stress  is  relatively  insensitive  to  the  core  thickness.    However, 
"o?      1 :      -est  to  the  contact  surface  tncreases,  the  shear  stress  transtuon 
mlg  li  «  the  tnterfa'c,  w,th  the  core.    The  resuit  is  a  reduction  ,n  the  s  ea 
I;:  atthe  interface  for  lavups  w.th  th.c.er  extenor  lamma.   it  ts  also  noteworthy  that 
lhe  maxtmum  shear  stress  occurs  at  cross  sect.on  B  vice  cross  sect.on  A. 

F.sures  36  37  and  38  are  graphs  for  the  transverse  normal  stress.  .,  .  for  the  same 
J—  layups.    These  graphs  show  an  increase  in  ..  as  the  th     nes s  of* 
,     o    lavers  increases     The  increased  th.ckness  of  these  layers  produces  a  stronger 

stresses  increase.    At  cross  sect.on  C,  some  tens.le  transverse  stress  ,s  euden,    A.  P« 
Is,  stated.  tens„e  transverse  stress  is  a  potentia,  source  of  delam.nat.on.    As 


56 


STRESS-Y:   ZOOM  OF  CONTfCT  ZONE 
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Figure  32.      Stress  contour:  oy  for  clamped  beam  model  (In  region  of  contact) 

0(8)-ISO(6)-0(2)  case.  The  resulting  smaller  contact  zone  lead  to  higher  contact  stresses. 
As  seen  in  the  symmetric  cases,  an  increase  in  the  thickness  of  the  layer  closest  to  the 
contact  zone  resulted  in  an  increase  in  the  magnitude  of  the  maximum  tensile  transverse 
stress,  seen  at  cross  section  C.  However,  the  stress  at  the  laminate  interface  was  always 
compressive. 

With  regard  to  bending  stresses  for  these  layups,  Figures  47  and  48  clearly  show  the 
nonlinear  behavior  as  the  contact  layer  thickness  increases.  In  addition,  the  thickness 
of  the  exterior  layer  opposite  to  the  contact  surface  shows  similar  results  as  the  sym- 
metric cases.  As  the  thickness  of  this  layer  decreases,  the  beam  is  more  susceptible  to  a 
bending  crack  that  propagates  into  the  interface. 
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Figure  33.      Stress  distribution:  zxy  for  0(3)-ISO(10)-0(3)  laminate 
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Figure  34. 


Stress  distribution:  t     for  0(4)-ISO(8)-0(4)  laminate 
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Figure  35.      Stress  distribution:  r    for  0(5)-ISO(6)-0(5)  laminate 
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Figure  36.      Stress  distribution:  ay  for  0(3)-ISO(10)-0(3)  laminate 
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Figure  37.      Stress  distribution:  <ry  for  0(4)-ISO(8)-0(4)  laminate 
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Figure  38. 


Stress  distribution:  <t   for  0(5)-ISO(6)-0(5)  laminate 
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Figure  39.      Stress  distribution:  <jx  for  0(3)-ISO(10)-0(3)  laminate 
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Figure  40.      Stress  distribution:  ax  for  0(4)-ISO(8)-0(4)  laminate 
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Figure  41.      Stress  distribution:  <r,  for  0(5)-ISO(6)-0(5)  laminate 
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Figure  42.      Stress  distribution:  ax  for  0(5)-ISO(6)-0(5)  laminate  at  cross  section  C 
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Figure  43.      Stress  distribution:  t,,  for  0(3)-ISO(6)-0(7)  laminate 
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Figure  44.      Stress  distribution:  xxy  for  0(8)-ISO(6)-0(2)  laminate 
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Figure  45.      Stress  distribution:  ay  for  0(3)-ISO(6)-0(7)  laminate 


72 


□ 

w 1 

o 

■ •— 1 

•   1 

n 

o 

*f 

r 

o 

T 

O 

i  n;rwn 

T 

O 

o-  c 

«-  C 

•  -  c 

ROSS  SECTION 
ROSS  SECTION 
ROSS  SECTION 

R 
B 
C 

1 

o 

ri 

0.0 


0.2  0.4  3.8  0.8 

NORMAL I ZED  DEPTH  BELOW  CONTRCT  SURFACE 


L.O 


Figure  46.      Stress  distribution:  <ry  for  0(8)-ISO(6)-0(2)  laminate 
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Figure  47.      Stress  distribution:  <rx  for  0(3)-ISO(6)-0(7)  laminate 
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Figure  48.      Stress  distribution:  ax  for  0(8)-ISO(6)-0(2)  laminate 
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V.     CONCLUSIONS  AND  RECOMMENDATIONS 

This  study  has  developed  two  methods  for  approximating  contact  stresses  using  the 
augmented  Lagrange  multiplier  method.  As  illustrated  in  Part  A  of  Chapter  IV,  these 
methods  accurately  approximate  the  stresses  that  result  from  contact  between  a  cylinder 
and  plane  surface.  This  study  has  also  illustrated  how  this  approach  can  be  applied  to 
understand  the  behavior  of  an  actual  contact  problem  by  examining  the  response  of  a 
composite  plate  to  low  velocity  impact. 

A.     RAYLEIGH-RITZ  APPROACH 

In  the  process  of  developing  these  methods,  a  number  of  comments  can  be  made 
regarding  the  application  of  the  Rayleigh-Ritz  method  to  solving  contact  stress  prob- 
lems. 

1.  The  selection  of  the  trial  function  is  an  extremely  challenging  process.  If  it  is  de- 
sired to  determine  the  deformation  in  a  contact  problem,  the  proper  stress  field 
must  be  first  satisfied.  Because  of  this,  the  selection  of  possible  trial  functions  is 
limited.  For  example,  when  selecting  a  trial  function  for  the  vertical  deformation 
of  the  foundation  of  Figure  1.  a  suitable  selection  is  given  by  the  following 
equation: 


v{xy)  =AX)  cos(  -fff  ) 


This  equation  exhibits  the  favorable  characteristics  of  maximum  deformation  at  the 
contact  surface  and  diminishing  deformation  as  the  distance  from  the  contact  sur- 
face increases.  If  contact  stresses  are  to  be  modeled,  this  trial  function  is  inappro- 
priate.   Calculation  of  ty  is  as  follows: 


-  =  -/Wl7Tsin(ltf  ) 


~  dy 


This  function  exhibits  zero  strain  at  the  point  of  contact  increasing  to  maximum 
strain  at  the  lower  boundary. 

2.  Since  the  selection  of  trial  functions  is  difficult,  the  task  is  further  impeded  by 
complicated  geometries.  Furthermore,  selection  of  a  trial  function  necessitates  that 
some  knowledge  of  the  deformation  field  exists.  Without  a  sensible  selection  of 
trial  functions,  an  accurate  approximation  is  unlikely. 

3.  This  method  assumes  the  trial  function  in  the  form  of  an  infinite  series.  Solution 
accuracy  theoretically  should  improve  with  an  increased  number  of  terms.  How- 
ever, precautions  must  be  taken  to  ensure  the  solution  is  numerically  stable  as  the 
number  of  terms  increases.  Since  the  strain  energy  calculations  require  integration, 
there  are  choices  of  trial  functions  that  will  increase  without  bound  as  the  number 
of  terms  is  increased.  This  problem  can  be  controlled  by  normalizing  dimensions 
or  limiting  the  choice  of  trial  functions. 
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4.  An  increase  in  accuracy  was  observed  as  the  number  of  constraints  was  increased 
and  the  distance  between  consecutive  constraints  was  decreased.  It  is  believed  that 
the  improved  accuracy  results  from  a  better  definition  of  the  contact  surface. 

B.  FINITE  ELEMENT  APPROACH 

The  results  in  Chapter  III  illustrated  that  this  approach  of  applying  the  finite  ele- 
ment method  to  contact  stress  analysis  is  effective.  A  number  of  comments  can  be  made 
regarding  this  approach  to  problem  solving. 

1.  As  illustrated  in  the  results,  this  method  accurately  approximated  the  isotropic 
roller  bearing  problem.  However,  some  difficulties  were  encountered  during  the 
modeling  of  the  multi-ply  composite.  In  this  model,  smooth  trends  of  decreasing 
deformations  were  often  interrupted  by  spurious  deformations  or  groups  of  defor- 
mations. These  interruptions  occurred  within  layers  of  significantly  reduced 
stiffness.  It  is  believed  that  the  optimization  routine  had  difficulty  approximating 
the  deformations  through  these  layers  because  of  their  very  small  contribution  to 
strain  energy.  As  stated  in  the  results,  the  contact  boundary  conditions  were  ob- 
tained from  the  optimization  program  and  applied  to  a  direct  finite  element  pro- 
gram to  solve  the  problem.  The  above  difficulty  is  recognized  as  a  limitation  of  this 
approach. 

2.  This  approach  is  much  more  flexible  for  complicated  geometries  than  the 
Rayleigh-Ritz  approach.  In  addition,  detailed  knowledge  of  the  deformation  field 
is  not  needed  as  required  by  the  Rayleigh-Ritz  approach. 

3.  The  application  of  static  condensation  is  crucial  to  the  successful  implementation 
of  this  method.  Every  effort  should  be  made  to  reduce  the  number  of  design  vari- 
ables to  improve  optimization  efficiency. 

C.  COMMENTS  ON  OPTIMIZATION 

A  number  of  observations  were  made  regarding  the  general  use  of  the  Automated 
Design  Synthesis  System  and  the  specific  usage  of  the  augmented  Lagrange  multiplier 
method. 

1.  The  global  optimum  was  more  likely  to  be  determined  when  the  objective  function 
was  normalized. 

2.  For  both  the  Rayleigh-Ritz  approach  and  the  finite  element  approach,  the  initial 
choice  of  design  variables  had  a  significant  affect  on  the  possibility  of  obtaining  the 
global  optimum.  For  the  Rayleigh-Ritz  method,  initial  selections  of  design  vari- 
ables can  result  in  largely  dissimilar  values  of  strain  energy  and  external  work,  the 
two  components  of  the  objective  function.  It  was  determined  that  convergence  was 
more  likely  when  optimization  commenced  with  these  two  terms  on  the  same  order 
of  magnitude.  With  regard  to  the  finite  element  approach,  sensible  choices  of  the 
initial  design  variable  vector  was  necessary  for  convergence  to  the  global  optimum. 
This  was  accomplished  by  intuitive  selection  of  design  variables  to  model  the  likely 
deformation. 

3.  Solution  accuracy  can  be  improved  by  scaling  constraint  equations.  It  has  been 
stated  that  in  some  circumstances,  some  constraints  change  more  rapidly  than 
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others  and  can  influence  the  solution  excessively  while  others  have  little  influence 
[Ref.  6:  p.  136]. 

4.  With  regard  to  the  usage  of  the  augmented  Lagrange  multiplier  method,  it  was 
frequently  necessary  to  tune'  the  optimization  algorithm  to  a  specific  problem. 
This  was  done  by  varying  the  initial  penalty  term  p  and  the  initial  Lagrange  multi- 
plier term  /  .  As  stated  by  Vanderplaats,  commencing  optimization  with  a  small 
value  of  p  should  theoretically  suffice  for  most  problems  [Ref.  6:  pp.  137-138]. 
However,  it  was  frequently  necessary  to  select  an  initial  value  for  p  due  to  conver- 
gence to  unrealistic  solutions.  Similarly  ,  an  initial  choice  of  the  Lagrange  multi- 
plier term  can  effect  the  solution.  Commencement  with  a  small  value  is  again 
recommended.  [Ref.  6:  p.  141].  This  need  to  tune'  the  problem  is  a  significant 
drawback  to  using  this  optimization  method.  The  ideal  way  to  overcome  this  lim- 
itation is  to  first  tune  the  optimization  routine  using  a  known  solution.  With  this 
accomplished,  this  approach  can  be  used  for  meaningful  data  collection. 

D.  SANDWICH  COMPOSITE  MATERIAL  STUDY 

The  behavior  of  sandwich  composite  materials  to  low  velocity  impact  loading  was 
successfully  investigated  by  the  application  of  the  finite  element  approach.  A  number 
of  observations  can  be  made  from  examining  the  results. 

1.  The  maximum  shear  stress  is  relatively  insensitive  to  layer  thicknesses.  However, 
as  the  thickness  of  the  contact  layer  increases,  a  reduction  of  the  interface  shear 
stress  is  observed. 

2.  Tensile  transverse  normal  stresses  exist  at  some  cross  sections  away  from  the  con- 
tact zone.  However,  this  stress  is  always  compressive  at  the  interface.  Compressive 
transverse  stresses  increase  in  beams  with  smaller  cores  due  to  reduced  deflection 
and  contact  zone  size. 

3.  As  the  thickness  of  the  layer  closest  to  the  contact  zone  increases,  a  nonlinear  dis- 
tribution of  bending  stress  within  this  layer  intensifies.  This  phenomenon  is  local- 
ized to  the  region  of  contact. 

4.  As  the  thickness  of  the  layer  opposite  to  the  contact  zone  increases,  bending  crack 
propagation  toward  the  core  is  less  likely  due  to  increased  compressive  bending 
stresses  within  the  layer. 

E.  RECOMMENDATIONS  FOR  FURTHER  STUDY 

The  methods  developed  in  this  study  offer  a  basis  from  which  additional  research 
can  grow.  A  reasonable  direction  is  the  relaxation  of  some  of  the  assumptions  made  in 
Chapter  II  Part  B.  For  example,  relaxation  of  the  rigid  roller  assumption  and  the 
frictionless  surface  assumption  would  provide  challenging  research.  Models  with  com- 
plex geometry  could  be  created.  For  example,  a  model  of  a  pin  loaded  bolt  connection 
could  be  created  with  rigid  or  non-rigid  pins.  Implementation  of  these  changes  would 
provide  a  versatile  and  highly  applicable  model  for  contact  stress  analysis. 
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